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1 .0 SUMMARY 


Described herein is a computer program for the analysis of the sub- 
sonic, swirling, reacting turbulent flow in an axisyimnetr ic , bluff-body 
research combustor. The program features an improved finite-difference 
procedure designed to reduce the effects of numerical diffusion and a new 
algorithm for predicting the pressure distribution within the combustor. 

A research version of the computer program described in the report was 
supplied to United Technologies Research Center by Professor A. D. Gosman 
and his students, R. Benodeker and R. I. Issa, of Imperial College, London. 
The Imperial College staff also supplied much of the program documentation 
contained in this report. 

This report presents a description of the mathematical model for flow 
within an axisymmetric bluff-body combustor, the development of the finite- 
difference procedure used to represent the system of equations, an outline 
of the algorithm for determining the static pressure distribution within 
the combustor, a description of the computer program including its input 
format, and the results for representative test cases. 

This report constitutes the final report for Task VII of NASA Lewis 
Research Center Contract NAS3-22771. Funding for this effort was provided 
by the Air Force Aero Propulsion Laboratory under contract FY1455-82-N0633 . 
The APAPL program manager was Dr. W. M. Roquemore. 

The NASA project monitor is Dr. C. J. Marek. 
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2.0 INTRODUCTION 


This report presents a description of the mathematical basis for and 
operation of a new version of the TEACH computer program developed 
originally by workers at Imperial Colleges London. The program is intended 
for use in the analysis of the subsonics swirling, reacting, turbulent flow 
in a cylindrical bluff-body research combustor (Fig. 1). This combustor 
has been developed by V. M. Rocquemore and colleagues of AFAPL as a 
research tool for gas turbine combustor modelling and for diagnostic 
instrumentation development (Ref. 1). The flowfield of interest is that 
which is developed in the open-ended, cylindrical chamber downstream of the 
base of the bluff-body which houses the fuel injector (Fig. 1). Such flows 
can be described by the fully-elliptic , steady-state equations of motion. 
The present version of TEACH solves these equations using an improved 
finite-difference procedure—the Bounded Skew-Upwind Differencing (BSUD) 
method. The axisymmetric pressure distribution is estimated by means of a 
new algorithm — The Pressure-Implicit Split Operation (PISO) predictor- 
corrector technique. 

Sections of the computer program have been reorganized by UTRC in an 
attempt to minimize the changes necessary to install the program on diffe- 
rent computer systems. The input format has been extensively revised to 
permit more flexibility in setting up and running cases. In the program 
that was supplied to UTRC, Imperial College included a number of additional 
options. These include switches that permit (1) the analysis of two- 
dimensional (planar) flows and (2) the selection of the more conventional 
hybrid differencing scheme. Verification of the operational status of 
theje optional features was beyond the scope of the present contract and 
would have required additional documentation from Imperial College. These 
options have been retained in the present version of TEACH for the 
convenience of the user. 

In Section 3, the system of equations describing the flow within the 
research combustor is presented. In Section A, the finite-difference forms 
of these equations, the special numerical treatment for certain source 
terms, and the PISO algorithm are described. In Section 5, the input 
format, source language listings, and other documentation for the computer 
program are presented. In Section 6, three representative sample cases are 
used to illustrate the results obtained with the computer program. 

References are listed in Section 7 and the nomenclature is presented 
in Section 8. 
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3.0 MATHEMATICAL MODEL 


Id this section } the governing equations describing the subsonic, 
Bvirling, reacting turbulent flow in an azisymmetric research combustor are 
presented. The tvo*equation turbulence model is then described. Finally, 
a combustion model for the turbulent reaction of propane vith air is 
outlined. 


3.1 Governing Equations 

The governing conservation equations for the mean motion of a two- 
dimensional (axisymmetric) , steady-state turbulent flow are presented below 
using cylindrical coordinates. Closure of this set of equations is 
provided by using the turbulent kinetic energy k and energy dissipation 
rate t which are derived from additional transport equations (see Sect. 3.2 
and Launder and Spalding (Ref. 2)) to obtain an effective eddy viscosity 
such that the unknown turbulence diffusional fluxes can be expressed as the 
product of the eddy viscosity and the gradient of the appropriate dependent 
variable; i.e., from the "gradient transport hypothesis" (Hinze, Ref. 3). 

The governing equations are: 

Continuity 


r~ (r:u) — (rrv) * 0 

ox dr 


( 3.1 .1 
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where represents such scalars as temperature, mass fraction, turbulent 
kinetic energy, etc. The quantities v and T are the effective 
exchange coefficients and are the sum of both tfie^laminar and turbulent 
transport coefficients 


■'eff 


‘"t* ^eff 


(3.1.6) 


Hi ^t 

-a. 4 

♦t 

The sources S, in the momentum equations represent additional 

terms associated with non-uniform viscositys Their influence is generally 
small except where changes in fluid properties have considerable effects. 
These terms are given by: 
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The governing equations (3.1.1) to (3.1.5) can be represented by the 
following general equation: 


where it is now assumed that the axial and radial pressure gradients have 
been incorporated into and S^, respectively. 


3.2 Turbulence Model 


pie two equation (k-c) turbulence model, developed by Launder and 
Spalding (Ref. 2), is used in the computer program. In the model, the 
turbulent viscosity is determined from the time-mean values of the kinetic 
energy of turbulence (k) and the volumetric turbulent kinetic energy 
dissipation rate ( c) by the relationship: 


V 


t 



(3.2.1) 


where is a constant. The quantities k and c are determined from 
transport-equations of the same form as Eq. (3.1.10) where for 
the k equation: 
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(3.2.2) 


« - k, S • V G “ C. r^k^/L 
^ t t 

and for the t equation 

« - c , S{ * Cl G p k - C2tE^/k 


where 




(3.2.3) 


( 3 . 2 . 4 ) 


The quantity G represents the production of turbulent kinetic energy by 
the mean motion. Turbulent kinetic energy is dissipated at a rate 
proportional to c which, in turn, is also generated and dissipated at rates 
that depend on k and z as can be seen from £q. (3.2.3). Thus, the 
equations for k and c are strongly coupled. In the vicinity of the walls, 
the determination of both k and t requires special consideration as noted 
in Section 4.3. The constants appearing in the turbulence model are listed 
in Table 1. 


3.3 Combustion Model 


The model proposed by Magnussen and Hjertager (Ref. 4) is used in the 
computer program. This model relates the rate of combustion of the fuel 
(restricted to propane in the present analysis) to the rate of dissipation 
of turbulent eddies and expresses the rate of reaction in terms of the mean 
concentration of the reacting species, the turbulent kinetic energy, and 
the rail of dissipation of turbulent kinetic energy. 


For propane, the stoichiometric relationship is: 

C 3 Hg 4 502 - 3 C 02 ♦ 4 H 2 O + 5.006x10^ J/kg (3.3.1) 

In the combustion model, the species concentrations are determined from the 
governing partial differential equations for the fuel mass fraction and the 
mixture fraction defined below. The enthalpy is defined as: 


h 


“f H 


Hd 4 C 


p,m 


(3.3.2) 


The mean specific heat is obtained from 


C 

p,m 


? *i <=pi 


(3.3.3) 


vith the summation over all species. The specific heats of individual 
species are obtained using polynomials of the form: 

Cp.i “ ^ * ^4iT^ ♦ ^3.3.4) 

where R is the universal gas constant, 8314.3 J/kg-molR). 


The constants of the polynomials are obtained from standard thermochemical 
tables and are listed in Table 2. 
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The mixture fraction (f), enthalpy (h) and fuel mass fraction (m,) are 
governed by equations identical in form to Eq. (3.1.10). For both the mix- 
ture fraction and enthalpy, the source term is zero. For the fuel mass 
fraction, the source term is calculated in accordance with the Magnussen- 
Hjertager model: 


- c f «in [ Amf, A “o2, AB ^ 1 

I ~T i*l\ (3.3.5) 

for which the constants are A ■ 4, B ■ 0.5; for propane, the stoichiometric 
oxygen to fuel mass ratio is 3.635. The rate of combustion is then, 


R£ • ~Sm£ 


(3.3.6) 


The fuel mass fraction is defined by the sum of the unburned and burned 
fuel : 


f 


1 

“f 3 "C02 


(3.3.7) 


where, from (Eq. 3.3.1) and the fact that propane and carbon dioxide have 
nearly identical molecular weights, the mass fraction of burned fuel equals 
one-third of the mass fraction of CO . Since both the mixture fraction and 
fuel mass fraction are determined by the solution of the respective 

transport equation, Eq. (3.3.7) can be used to determine the mass fraction 
for CO^. 


mc 02 * ^ (f-mf) 


(3.3.8) 


Similarly, for the remaining species: 

®H20 “ 0-545 mco2 

- 1.635 (f-mf) 


(3.3.9) 


®pr • “k :02 + mH 20 * ^-635 (f-m^) - (i+l)(f-ni^) (3.3.10) 


BN2 - 0.767 (1-f) 

mo2 - 1 - (0^2 ♦ Of ♦ mco2 ♦ mjj20^ 
■ 0.233 - 3.868f ♦ 3.635 mf 

« 0.233 - (i+0.233)f ♦ imf 
The average density is given by 

^ “ RT r 

i Mf 

and the temperature is calculated from Eq. (3.3.2). 


(3.3.11 ) 


(3.3.12) 


(3.3.13) 
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4.0 FINITE-DIFFERENCE APPROXIMATIONS 


In this section, the finite-difference approximations used in the 
computer program are discussed. In Section 4.1, the Bounded Skew-Upvind 
Differencing (BSUD) scheme is presented. In Section 4.2, the special 
treatment given to certain source terms in the equations of motion are 
discussed briefly. In Section 4.3, the boundary condition formulations are 
outlined. Finally, in Section 4.4, the predictor-corrector algorithm for 
determining the static pressure distribution is presented. 


4.1 The Bounded Skev-Upvind Differencing Scheme 

In this section, the Bounded Skev-Upvind Differencing (BSUD) scheme is 
described. First, a brief reviev of the flux form of the equations of 
motion is presented. Second, a detailed description of the finite- 
difference form of the flux contribution to a representative face of a 
typical control volume is given; the derivation of the flux contributions 
to the other faces is then outlined. Third, the resulting coefficients 
for the finite-difference equations representing the total flux (and 
sources) are presented. Finally, the bounding scheme for the coefficients 
is detailed. 


4.1.1 Flux Form of the Equations of Motion 

Further generalization of £q. (3.1.10) leads to the conclusion that 
the equations of motion for both laminar flow and (time-averaged) 
turbulent flow can be vritten in similar fashion for all of the dependent 
variable » : 


^ (puO 

0X 


3r 



7 ^ 


3r ' ^ 3r ' 


( 4 . 1 . 1 . 1 ; 


vhere f « 0 for two-dimensional (planar) flow and 6 » 1 for axisymmetric 
flow. The variable C represents any of the dependent variables (e.g., the 
velocity components u, v, w, mixture fraction, turbulent kinetic energy 
and turbulent energy dissipation rate, etc.). The exchange coefficient, 

(f # represents the sum of both laminar and turbulent contributions and is 
interpreted as the eftective viscosity for 4> - u, v, w, the effective 
diffusivity for < - mixture fraction, etc. S- is a generalixed source 
term. ^ 


Eq. (4. 1.1.1) is integrated over a control volume appropriate for each 
dependent variable C and, after some manipulation, the finite-difference 
equivalent form of Eq. (4. 1.1.1) is obtained. The control volumes are 
defined using an orthogonal grid formed by the intersection of co-ordinate 
lines in the axial and radial co-ordinate directions. The intersection of 
8^^i^ lines. Fig. 2, form the grid nodes at vhich all flow properties 
except the axial (u) and radial (v) velocities are calculated; i.e., all 
scalars and the tangential velocity component (w) for swirling flows. The 
axial velocity is calculated using a second grid with grid lines such that 
the grid nodes of the scalars in the axial direction are midway between the 
velocity grid nodes. Fig. 3. The axial velocity and scalar grid 
nodes are co-incident in radial position. The radial velocity is 
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calculated using a third grid such that the grid nodes of the scalars in 
the radial direction are midvay between the radial velocity grid nodes » 

Fig. 3. The radial velocity and scalar grid nodes are co-incident in axial 
position. Directions in the grids are identified as northi south, east and 
vest. It should be noted that it is the axial grid line locations for the 
axial velocity and the radial grid line locations for the radial velocity 
which are input to the present version of the computer program. 

Control volumes for each scalar and the tangential velocity are 
defined such that (1) the east and west faces of the control volume are co- 
incident with the axial velocity axial grid line locations and (2) the 
north and south faces are co-incident with the radial velocity radial grid 
line locations as shown in Fig. 4. The control volumes for the u velocity 
are shifted relative to the scalar control volumes such that the east and 
vest faces are co-incident with the scalar grid node axial locations; 
similarly, the v velocity control volumes are shifted such that the north 
and south faces are co-incident with the scalar grid node radial positions. 


The finite-difference form of Eq. (4. 1.1.1) is derived by integrating 
this equation over the appropriate control volume. In performing the inte- 
gration over the control volume for each term in £q. (4.1. 1.1), the mean- 
value theorem is employed and the source term is linearized in the vicinity 
of the center of the control volume (point P). After some manipulation, 
the finite-difference form of Eq. (4.1 .1.1) is obtained 

'l*. - Vv * S'n - V. ■ ■ *p’ ■ ®W “p - *w’ 

(4. 1.1. 2) 


♦ “l. % - *p’ • ®S “p • ‘s’ * * Vp’ 

ive c 

(ou). 


where C_, C , etc. are "convective coefficients" as defined below 

JL wt 


a 

e e 


H • *n 

Cv ■ <Pu>v *w 


(4.1.1 .3) 


Cg • (pv), «, 

and a « a , a ■ a are the areas of the faces of the control volume. The 
"diffusion* colfficients are given by 



and ('*) is the distance between points P and E, etc. (e.g.< see Fig. 4). 
e 


4-2 



It is important to note that £q* (A.1.1.2) applies to all of the 
dependent variables although the appropriate grid must be used in each case 
to define the geometric parameters used in the calculation. Also, Eq. 

(4. 1.1 #2) applies to both the hybrid and bounded skew-upwind differencing 
procedures used in the computer program since each scheme is simply an 
alternative method for defining (interpolating for) the fluxes at the faces 
of the control volume (e.g., C * ^ * However, the diffusion terms 

are always represented by central differences. 

It is convenient to define a total flux for each face of the control 
volume as the sum of a convective flux and a diffusive flux such that 


r, - r, ♦ r„ - r. - $„ ♦ <4. 1 . 1 . 5 ) 

where 

Fw “ 

(A.l .1 .6 ^ 

^ “ Vn ‘ ^ 

F, . Cgtf, - Dg Up - <5) 


In the following section, the skew upwind differencing procedure will be 
used to calculate the values of the dependent variables at the faces of the 
control volume. As a result, Eq. (4.1 .1.5) will include not only the 
values of at the •"normal", or main, grid node locations (E, W, K, S and 
P) but also at the corner locations (NE, SE, NW, and SV). The finite- 
difference form for Eq. (4. 1.1. 5) will be shown to be 


Vp * '^e'e ^s's " * ^se'se 

(4.1 .1 . 7 ) 


4.1.2 - Calculation of the Fluxes 

Recall that the equation of motion, Eq. (4. 1.1.1), can be written in 
terms of fluxes to each face of the control volume, Eq. (4. 1.1. 5). In this 
section, a procedure will be described to calculate fluxes, F , F , F^, and 
F . The derivation of F , the flux to the west face, of a typicaY scalar 
control volume is given in detail. The derivation of the fluxes for the 
other faces or the fluxes for the u and v velocity components are outlined. 

Consider the control volume shown in Fig. 5. For this case, it has 
been assumed that the velocity vector is oriented as shown; i.e., the u and 
V components are non-negative. The flux to the west face is given by: 

Fy • “ Dy (tp ” 4y) (4. 1.2.1) 
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For central-differencing (CD), the value of the dependent variable at the 
west face, is given by linear interpolation between ♦ at the W and P 
grid nodes : 


♦w • (I - 


(4. 1.2. 2) 


where the interpolating factor is 


_ Xw - Xw 

" V - V. 

*p *w 


(4. 1.2. 3) 


The central difference form of the flux at the west face is then; 


^*CD 

»W 


1'^, <>-“.>*» 1 >> 


(4. 1.2. 4) 


where the Peclet number at this face is given by: 

**«w " (4. 1.2. 5) 

The upwind difference (DD) form for the flux at the west face is obtained 
by setting to rero and neglecting the diffusion term: 

Dy " ^*w (4. 1.2. 6) 

Equation (4. 1.2.1} is also the starting point for the skew upwind 
differencing (SUD) scheme. The value of at the west face of the control 
volume, e^, is determined by extrapolating the velocity vector upstream to 
the point w' which lies along the grid line connecting the west and 
southwest nodes (see Fig. 5) to give: 

- (1 - kv> #y ♦ kw ♦sw (4. 1.2. 7) 

where the skew interpolation factor is the ratio of the vertical distance 
between w' and SW to the vertical distance between W and SW: 

w .1 h M 

^ 2 Ay (4. 1.2. 8) 

For very large flow angles (skewing) relative to the co-ordinate 
directions, k^ will exceed unity and i^w' will be defined in terms of C at 
the SW and S nodes; however, it is known that this approach can yield 
negative coefficients at the corner nodes (HW, SW, BE, SE) which can in 
turn produce oscillations in the solution. To assure that the coefficients 
for the corner nodes are non-negative, then: 
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The use of absolute value in Eq. (4. 1.2. 9) permits this equation to be used 
to define k for all velocity components at the vest face. 


In bis original development of the skev upwind differencing approxima- 
tion, Raithby (Ref. 5) assumed that i ■ . Thus, 

w w 

” '^p “ (4.1.2.10) 

It is desirable to use the central-difference procedure for small values of 
the grid Peclet number and the skev upwind differencing method for large 
values of the grid Peclet number# It is also desirable that these two 
formulations produce a continuous transition at the transition Peclet 
number which in the present case is: 





At the transition Peclet number, the central difference result (Eq* 

4. 1.2. 4) is: 

_ 5 V 

I>W “ (4.1.2.12) 


while the skew upwind differencing method (Eq. 4.1.2.10) yields: 




— - - ♦w) 

w 


( 4 . 1 . 2 . 13 ) 


Noting the definition given by Eq. (4. 1.2. 7), it is clear that these two 
results are not equal. 


The fluxes at the transition Peclet number can be made equal by noting 

(contrary to the assumption made by Raithby) that < and c , are related 
by: w w 


♦w' 



As ♦ 


ao that Eq. (4.1.2.10) becomes 


(4.1 .2.14) 


^*SUD 

I>W 





As - (4p - 


(4.1.2.15) 


Writing the central difference result in terms of the flux definition, Eq. 
(4. 1.2.1) then : 



(4.1.2.16) 


^ - <♦» - *w> 

Clearly,, these two fluxes will be equal at the transition Peclet number if 
a correction Pe*0((^/oS) As is added to the akev upvind differencing flux, 
Eq. 4.1.2.10, to obtainY 


^*>SUD 


Pe^4„. 



4. - (ip - ♦v,) 


(4.1 .2.17) 


The derivative (°'/^s)v can be computed by 

i± . *- ~ ♦»' 

*» As (4.1.2.18) 

Then, using Eq. (4. 1.2. 2), (4. 1.2. 7) and (4.1.2.11), Eq. (4.1.2.17) 
becomes 


^*SUD 




tPe^ - \ <*W ■ ♦sw’ 


(4.1 .2.19) 


At the transition Peclet number, the fluxes calculated by central 
differencing (Eq. 4.1.2.12) and skew upwind differencing (Eq. 4.1.2.19) 
are equal. 

It will be recalled that the above result for the skew upvind 
differencing flux at the vest face of the control volume was derived for 
non-negative values of the axial and radial velocities. Similar results 
can be derived for other combinations of u and v by consistent application 
of the process leading to Eq. (4.1.2.19). The result is a general 
expression for the flux as calculated using skew upvind differencing: 


^**SUD 


■ ♦W ♦ ♦pi 

- ^> ( < { 

» (1 - - o’ }) 


(4.1.2.20) 


%x ▼ 

The parameters, o and o , are switches that indicate the direction of the 
components of the'^local Ylov velocities 


• i * to) 


(4.1.2.21) 
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(4.1.2.22) 


Each of these parameters has a value of unity if the velocity component is 
positive (or, by convention, non-negative) and aero if it is negative. The 
transition Feclet number is now given by the general result. 


P 


«v 


"TT 

0 


1 

(I - a*) 


(4.1.2.23) 


In the hybrid differencing procedure, the more accurate central 
differencing formulation (Eq. 4. 1.2. 4} is used when the Peclet number is 
less than the transition value while the leas accurate, but stable, upwind 
differencing result (the generalization of Eq. (4. 1.2. 6)) 


I>W 




(1 - o;;) 


(4.1.2.24) 


is used when the Peclet number is greater than the transition value. Ori- 
ginally it was believed that a similar hybrid procedure could be developed 
for skew upwind differencing with Eq. (4. 1.2. 4) used for Pe < Pe* and Eq. 
(4.1.2.20) used for Pe > Pe* . However, this approach proved to be 
unworkable since some of the coefficients derived from this hybrid 
formulation for use in Eq. (4. 1.2. 7) can be negative. As an alternative, a 
flux blending scheme is used in which a weighted average of the upwind 
differencing and the skew upwind differencing fluxes is used. The 
weighting (blending) factor, \, is chosen in such a way as to assure 
boundedness (i.c., all co-efficients of Eq. (4.1.1 .7) are non-negative). 
The bounded skew-upwind differencing (BSUD) flux is defined by 


" T*,) (4.1.2.25) 

vith the veighting factor restricted to the range, 0 ^ As a conse- 

quence of this definition: 

[ °w \ ♦sw - dww} (4.1.2.26) 

♦ (1 - o“) { dp - o;; - <i - oj) ♦„| ] 

Finally, a bounded skew hybrid differencing (BSHD) formulation can be 
defined as: 


^•BSUD " VwcD ♦ <1 - V ^nSUD (4.1.2.27) 


where Y is permitted to assime only two values: 't' ■ 1 for central 
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differencing (Pe < Pe*) and X ■ 0 for bounded skev-upvind differencing (Pe 
^ Pe* )• Eqs. (4.1.2.26) and (4.1.2.27) are the basic working 
relationships to determine the flux at the west face of the control 
volume. The contributions of the flux to each of the co-efficients in Eq. 
(4. 1.1. 7) can be immediately identified by using Eqs. (4.1.2.26) and 
(4.1.2.27) in Eq. (4. 1.1. 5). 

Equations analogous to Eq. (4.1.2.26) can be derived for the other 
three faces of the control volume in exactly the same manner as employed 
above. However, the results can be obtained by inspection as follows: 

^a^t 

The east face flux is obtained by translating the nodal subscripts 
eastward such that: 


West Subscript East Subscript 


w 

becomes 

P 

p 

becomes 

E 

sw 

becomes 

S 

NW 

becomes 

N 

s 

becomes 

SE 

N 

becomes 

NE 

Of course, the lower case 

subscript 

"w" becomes "e". 




The south face flux is obtained from the 

west flux Eq. 4.1.2.26 by 

rotating the nodal subscripts through 90 degrees in the counterclockwisi 

Vest Subscript 


East Subscript 

w 

becomes 

s 

sw 

becomes 

SE 

NW 

becomes 

SW 

S 

becomes 

E 

N 

becomes 

W 

o“ 

becomes 

0^ 
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1 U V 

X ^ becomes 1 - c 

becomes 1 - 

1 - becomes 

North 

The north face flux is derived by translating the south flux result 
northward : 


South Subscript Worth Subficript 


s 

becomes 

p 

p 

becomes 

N 

SE 

becomes 

£ 

sv: 

becomes 

V 

E 

becomes 

N£ 

W 

becomes 

NW 


The results of the manipulations are summarized in Tables 3 through 5. 

4.1.3 - Cal culation of the Coefficients for the 
Finite-Dif ference Form of the Equations of Motion 

The finite-difference form of the equations of motion, (e.g., Eq. 

(4.1. 1.7)) can be derived directly from the flux information presented in 
Tables 3 through 5 and the sign conventions determined from Eq. (4. 1.1. 5). 
The resulting expressions will contain the unknoira blending factor, \. The 
blending strategy requires that the terms in the equations for the coeffi- 
cients most responsible for producing negative co-efficients be isolated 
so that appropriate values for Y can be determined. Furthermore, the 

control volumes adjacent to the physical boundaries of 
the flow may require modification to incorporate the effect of the boundary 
conditions. Thus, to simplify manipulation and modification, some 
additional notation will be defined. 


Let the center of the control volume (point P) be located at the 1th 
axial position and Jth radial position. The flux contributions (the 
components of the total flux) to the east face are denoted as E1(I,J), 
E2(1,J), 0(1, J) and the flux contributions to the north face are denoted 

as N1(I,J), N2(I,J), N3(I,J). These flux contributions are defined as 
follows: 
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1 - o'* 

becomes 

1 - o'* 

c V 

becomes 

1 - o“ 

1 - 

becomes 

c« 


Hor.th. 

The north face flux is derived by translating the south flux result 
northward: 

South Subscript North Subscript 


s 

becomes 

P 

p 

becomes 

N 

S£ 

becomes 

£ 

SW 

becomes 

W 

E 

becomes 

N£ 

W 

becomes 

NV 


The results of the manipulations are summarized in Tables 3 through 3. 

4.1.3 - Calculation of the Coefficients for the 
Finite-Difference Form of the Equations of Motion 

The finite-difference form of the equations of motion, (e.g., Eq. 

(4. 1.1 .7)) can be derived directly from the flux information presented in 
Tables 3 through 5 and the sign conventions determined from Eq. (4. 1.1. 5). 
The resulting expressions will contain the unknown blending factor, Y. The 
blending strategy requires that the terms in the equations for the coeffi- 
cients most responsible for producing negative co-efficients be isolated 
so that appropriate values for can be determined. Furthermore, the 
coefficients for the control volumes adjacent to the physical boundaries of 
the flow may require modification to incorporate the effect of the boundary 
conditions. Thus, to simplify manipulation and modification, some 
additional notation will be defined. 

Let the center of the control volume (point P) be located at the Ith 
axial position and Jth radial position. The flux contributions (the 
components of the total flux) to the east face are denoted as E1(I,J), 
E2(I,J), E3(I,J) and the flux contributions to the north face are denoted 
as N1(I,J) N2(I,J), N3(I,J). These flux contributions are defined as 
follows : 
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^entr.al_ Differencing : 

E1(1,J) - De - o^C£ 

E2(I,J) • El(l.J) ♦ C£ 

E3(I,J) - 0 

Hl(l.J) - Dk - a„CK 

N2(I,J) - ♦ Cj; 

K3(1,J) ■ 0 

JBo^n^e^ Differencing 

EUl.J) • (c“ - 1 ) Cj 
E2(1,J) - El(I.J) ♦ Cj 
E3(I,J) « kgCCj- - Pe*I>E^ 

KKI.J) - (c^; - 1) c^, 

K2(I,J) - N1(I,J) ♦ Ci; 

N3(1,J) - k„(CK - p/^D^,) 

The use of central vs. bounded skew-upvind differencing is determined by 
the value of the Peclet number at each face. The parameters : ■ 

k^, ... are local values; the subscripts (1,J) have been omitted in^the 
interest of readability. The corresponding flux contributions at the vest 
face are given immeciately by E1(1-1,J), E2(1-1,J), E3{1-1,J) and the fi-x 
contributions at the south face are N1(I,J-1), N2(I,J-1), H3(I,J-1). 

The coefficients of the finite-difference form of the equations cf 
motion may then be defined in terms of these flux contributions. The 
results are presented in Table 6. The coefficient Ap can be shown to be 
equal to the siun of the other eight coefficients when use is made of the 
mass continuity restriction: 

^E “ “ Cg - 0 (4. 1.3. 5) 

The boundary conditions (Section 4.3) and the bounding scheme (Section 
4.1.4)) can be applied directly to the flux contributions so that the 
results shown in Table 6 are general. 

4.1.4 - The Bounding Scheme 

The calculation of the bounded skew-upvind differencing fluxes and, 
therefore, the determination of the coefficients to the finite-difference 
form of the equations of motion requires that the blending factor, >, be 
determined. The blending factor specifies the relative proportions of the 
flux computed using skew and upwind differencing. For example, at the west 
face of the typical control volume, Eq. (4.1.2.23) states: 


(4.1 .3 .1 ) 


(4.1 .3 .2, 


(4. 1.3 .3) 


(4.1 .3 .4} 
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^^BSUD 




( 4 . 1 . 4 . 1 ) 


The co-efficients including the local blending factor are listed in Table 

6 . 


It is possible to show that the corner co-efficients A. , A , 

unconditionally non-negativee For example, consider the 

coefficient A_,: 

sw 


^sw 


(I.J) . „„ ♦ r. ivD<3i,,j.i 


( 4 . 1 . 4 . 2 ) 


Both and ) g are restricted to the range 0 £ V £ 1. If both the axial 
and radial velocities at the vest face of the control volume are positive, 

then both ^ and at this location are unity; in any other case, either 
(or both) are zero. For positive u and v velocities, the flux contribution 
E3(I-1,J) at the west face is positive. Therefore, the west face flux con- 
tribution to Ag„is always nonnegative. The same reasoning when applied to 
the south face leads to a similar conclusion. Therefore, A is uncondi- 
tionally non-negative. The other corner coefficients are treated in the 
same manner . 


Now, consider the four main co-efficients (A,., A_, A_, A.,). Froc 
Table 6: T S Tv 


A^;(1,J) - E2(I-1,J) - y tc“E3], , , 

♦ Y, c“k3] 

n n 1,J-1 


( 4 . 1 . 4 . 3 ) 


- 'Tn lr»] 


V u^ 

c 

n n I,J 


From the definitions of the flux contributions, Eqs. (4. 1.3.1) through 
(4. 1.3 .4), it is evident that: 

(1) E2(I-1,J) is positive if the axial velocity is positive and it is 
zero otherwise; 

(2) 0 ^ and E3 at the west face (1—1, J) are poaitive if the axial velo- 
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city is positive but is zero if u is negative; 

(3) froiD the definitions of E2 and E3 and the fact that £ 1, it is 
therefore concluded that the first two terms in Eq. (4. 1.4. 3) must yield a 
Don*Degative result ; 

(4) the third term is zero unless u is positive and v is negative at 
the south face, in which case, it is negative; 

(5) the fourth term is non-positive because the term within the 
brackets is always non-negative. 

Therefore, it is possible that the co-efficient is negative. 

It is desirable that all of the co-efficients of the finite-difference 
form of the equations of motion (Eq. (4.1 .I,?)} be non-negative for in this 
case the value of the dependent variable at the node P is simply a 
weighted average of the values of c at the surrounding nodes exclusive of 
the (somewhat complicating) effects of local sources. A bounding scheme is 
a procedure to limit the values of the co-efficients of the finite- 
difference equations in such a manner as to produce this physically realis- 
tic result. Its principal computational advantage is to exclude under- ar.c 
overshoots of the solution during the iterative procedure; these oscilla- 
tions can produce severe numerical instability. 

The bounding procedure used herein is based upon the following 
sequence : 


(1) For convenience, the following quantities are defined: 


a - E2(I-1,J) - e3]^ ^ , 

e l,j-l 

(4.1 .4.^; 

= c“ N3]j J.J 

(4.1 .4.5) 

•n • f°n On *^3]^ ^ 

(4.1 .4.6) 

•o that 


*w - ! - 1-.S. - ■>„«„ 

(4. 1.4. 7) 

From the previous discussion, a > 0, i >0 and a >0 

s “ n 



(2) It is desired that A^ > 0. Therefore, the blending factors, 
and are the solution to the linear programming problem ® 


is>s ♦ Vn i « (4.1 .4.8) 
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subject to the constraints 


0 < > 


s 



(4. 1.4. 9) 




(3) The solution to this problem is given by: 


(4.1 .4.10) 


>. i.n 


* min [ Yg, i/a ] 


(4.1.4.11) 


max [0, (a - ) 


(4.1.4.12) 


if_a* > ks 


- min [ Y^, i/i„] 


(4.1.4.13) 


^5 * max [ 0, (a-a^-inl/Sgl 


(4.1 .4.14) 


The upper limits, and '• , are normally equal to unity or are the values 
determined for the faces of adjacent control volumes; e.g., the vest face 
of the control volume at (I, J) is the east face of the control volume at 
(I-l, J). 


The resulting set of blending factors are then used to compute all of 
the coefficients in Table 6. 


4.2 Source Terms 


In this section, the finite-difference approximation to representative 
source terms of the governing equations are described briefly. From Eqs. 
(4. 1.1. 2) and (4. 1.1. 7) it can be seen that the source term is linearized 
into a term independent (in an explicit sense) of the dependent variable 
at the grid node F and a term explicitly dependent on f : 


S 




S 

u 


♦ S ♦ 

P P 


(4.2.1) 


Then, using Eq. (4. 1.1.7), 


one "solves" for : 
P 




iAj<ti ♦ 

N.- Sp 


(4.2.2) 
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The determination of and S for each dependent variable may be 
illustrated by considering thf source term for the rate of reaction of 
fuel, Eq. (3.3.5): 


Sr.£ 


- r ^ min I ^‘f’ ^2> A®>"pr 
' i iTT 


(4.2.3) 


If either the fuel or oxidizer are the minimum of the three rates contained 
in the brackets, then: 


or 


Hence , 


Sil£ * - 



Sm£ * 0 



S 


u 


0 


(4.2.4) 


(4.2.5) 


(4.2.6) 


S 


P 


PcA 

k 


But if the third teru (the product term) is the minimum, then 


Sr.£ 


h - A Bripr 

k i + l 


(4.2.7) 


“ ■ AE (f-ir.£) (4,2.8) 

- - AE (f-2ir.£) - ^ ABr.£ 


Then, 


Sy - - AB (f-2r.£) 


(4.2.9) 


S 


P 



(4.2.10) 


It is seen that Eq. (4.2.9) depends explicitly on the fuel mass fraction; 
necessarily, is evaluated using the results from the previous iteration. 


Numerical stability in the determintion of is enhanced if the 
efficient S is non-positive. If S is positive.^then it can be seen 
the denominator in Eq. (4.2.2) may become negative (or aero) causing 
oscillations in the solution for . The linearisation of the source 


co- 

that 

term 
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for tfae fuel ness fraction guarantees that is non*-positive. 

Other source terms in the governing equations require similar 
treatment. For example, the source term (per unit volume) in the 
tangential (svirl velocity) momentum equation is: 

^ 2 ? 

- . - - _ (rvu.ff) 

which is written in finite-differences as: 


( 4 . 2.11 ) 


Sw 2 




I Weffn - rs«sWeff, 


( 4 . 2 . 12 ) 


It is maintained that the form of Eq. (4.2.12) is essential to preserving 
the equivalency between the partial differential equation for swirl 
velocity and the finite-difference approximation. Other forms such as: 




( 4 . 2 . 13 ) 


or 


<5V 


ro^r 


- (rvueff). 


( 4 . 2 . 14 ) 


do not preserve the desired equivalency because the quantities 
u -- and u ttust be interpolated at the desired locations (i.e., at 
ani^r ). xSe^source term, when multiplied by the cell volume (r .r.lx) 

u n » r ^ p 

becomes 


* - 2^xr^w^Ueff^ ^ 2lxrgVg^gf£^ 


( 4 . 2 . 15 ) 


or 


Sw • 24xrgWsU^j^^ ♦ 


(• 


2AxrnWniieff 




( 4 . 2 . 16 ) 


so that 


and 


Su “ 24xr,w^u^jj^ 

_ 2dxrnw„ueff„ 


( 4 . 2 . 17 ) 


(4.2.18) 


It must be noted that both w and w may depend on w and that the values 

s n p 
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of and v in Eqs. (4,2,17) and (4.2.18) are obtained from values at 

the previous iteration. 


4.3 Boundary Conditions 

The boundary conditions are applied to the flux contributions as 
defined by Eqs, (4.1 ,3.1 ) through (4.1 .3.4) so that the set of coefficients 
(e.g.. Table 6) for each dependent variable is the same throughout the 
computational domain. This consistency simplifies the application of both 
(1) the algorithm for solving the set of simultaneous equations for each 
variable and (2) the bounding procedure. 

The boundary conditions are applied to the control volumes in the 
vicinity of the physical boundaries for the flow field (see Fig. 1). Thus, 
the left-most control volumes include either the main flow or secondary 
flow inlets or the bluff body. The right— most control volumes include 
the combustor outlet. The top-most control volumes include a solid wall 
while the bottom-most control volumes include the axis of symmetry. 

^nl.et.s 

The flows in either the main or secondary inlet are known and 
specified. It is assumed that constant values of axial velocity, swirl 
velocity, turbulence intensity, fuel concentration, and temperature are 
known for each inlet. The radial velocity is assumed to be zero and the 
static pressure is assumed to be the same for both inlets. 

Axis_o^ Sym^t ry 


At the axis of symmetry, the normal gradient for each dependent 
variable, except for the radial velocity, is assumed to vanish; the radial 
velocity is assumed to be zero. By setting the co-efficients A^, and A tc 
zero for the radial gridline adjacent to the axis of symmetry, then the 
finite difference form of the normal gradient becomes: 



« 0 

Axis of symmetry 


(4.3.1) 


variables, except the radial velocity, the axis of symmetry is 
located midway between the J*1 and J«2 radial gridlines. For the converged 
solution, the boundary condition (Eq. (4.3.1)) produces the result: 

♦ l.l * <1,2 ) (4.3.2) 

where I is the «cial gridline index. For the radial velocity, the axis of 
aynmetry and the radial gridline for radial velocity are coincident and the 
boundary condition becones: 


V 


I.l 


(4.3.3) 
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At the outlet of the combustor > the axial gradient in all variables 
except axial velocity is assumed to vanish. By setting the co-efficients 
and A^ at the outlet to zero, the finite-difference approximation to the 
gradient becomes 

) -0 

/outlet (A. 3. 4) 

so that 

Vl.J " ^NI-1,J (4.3.5) 

where NI is the number of axial grid lines. For the axial velocity, a 
velocity correction is applied uniformly so that 

^NIU,J * ^NIl’-l,J ^INC (4.3.6) 

where NIU ■ NI-1 . The correction U is the velocity increment necessary 
to conserve mass flow at the outlet. 

SoUjd Wail 

If the resolution of the grid system could be made arbitrarily fine in 
the vicinity of a wall, then application of the appropriate boundary condi- 
tions at the wall would be straightforward. For example, the velocity com- 
ponents u, v, w vanish at a wall. However, for a realistic grid system, 
setting the velocity to zero for components of the velocity vector parallel 
to the wall (by adjusting the co-efficients of the finite-difference equa- 
tions) is improper because the co-efficients do not properly account for 
the effects of the wall; i.e., the shear is computed incorrectly. Instead, 
wall functions are used to compute contributions to the source terms in the 
governing equations. For a velocity component normal to the wall, the co- 
efficients may be adjusted so that this velocity component vanishes. (See 
the discussion for the axis of symmetry boundary condition.) 

Consider the case for which the axial velocity is parallel to a wall. 
The flow in the region near the wall is assumed to behave as one— dimension- 
al Couette flow. For turbulent flow, the shear stress is assumed to be 
constant in the Couette flow region except in the laminar sublayer 
immediately adjacent to the wall. Following Launder and Spalding (Ref. 2), 
the laminar sublayer extends from the wall to y-f « 11.63 where 

♦ I P C^ y \ 

y • I 1 ) (4.3.7a) 

P 

with the density, molecular viscosity, turbulent kinetic energy and grid 
distance from the wall evaluated at node P. If y+ is less than 11.63, then 
the shear stress is given by: 
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T « 

V 


(4.3.7b) 



P 


Otherwise the point P is within the turbulent flow region and the 
stress IS calculated using the logarithmic law of the wall so that 


shear 


c C. K u/£n 


(Ey*) J 


(4.3.8) 


Here K and E are constants 
to the source term for the 


given in Table 1 . The shear stress contributior. 
axial velocity momentum equation is 





(4.3 .9; 


where A- is the area of the face of the control volume adjacent to the 


The influence of the walls 
treated similarly. Of course, 
be interpreted appropriately, 
top wall must be resolved into 


on the v and w velocity components is 
the noraal distance y in Eq. (4.3.7) must 
For swirling flows, the shear stress on the 
axial and azimuthal components. 


For the turbulent kinetic energy equation, 
ponding to the wall node is set to zero and the 
is given by: 


the coefficient corres- 
source term (see Eq. 3.2 


2 .) 


S 


k 


V G 
t 



3 /: 

f k 



y"^ < 11.63 


(4.3 .1C; 


Sj^ ■ MiG 


3A 3/3 ^ 1 

C c k in (Ey^)/(Ky) 

I J p 


, y”^ 11.63 


(4.3.11 ) 


the generation term G to 

compute (uou/oy). 


The turbulent dissipation rate, c, is computed in the vicinity of a 
V.11 by a.aumiis th.t the rate, of turbulent kinetic energy generation and 


c 

p 


' 3/^. 3/2 ■ 

k 

Ky 

P 


(4.3.12) 
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(4.3.13) 


By defining 


the source for ^ as 


S 

c 


10 


30 


♦ 10 


30 


€ 

P 


then by Eq. (4. 1.1. 7), (t ■ c . 

P P 

For both the mixture fraction, f, and the fuel mass fraction, m , , it 
is assumed that the flux normal to the wall vanishes. Thus, for both f and 
nj this boundary condition is mathematically identical to the axis of 
symmetry boundary condition. 

The enthalpy at the wall is determined using Eq. (3.3.2) with the wall 
temperature specified and the fuel mass fraction calculated as above. 

Kail “ “fwall (4.3.14) 


4.4 Pressure Implicit Split Operator 


The finite-difference approximations to the governing equations are 
solved using the Pressure Implicit Split Opertor (PISO) developed at 
Imperial College by Issa (Ref. 6). The finite-difference equations are 
solved iteratively. The velocities are first calculated from the momentuc 
equations using a guessed pressure distribution! then the pressure distri- 
bution is adjusted so that the velocities satisfy the mass continuity 
equation and the cycle is repeated. 

The PISO algorithm involves splitting operations that couple the 
velocity and pressure variables. 

Let the superscripts *, ** and *** denote intermediate field values 
obtained during the operation of the algorithm which consists of a 
predictor step and two corrector steps. 

Pr^dic^or. Step - The pressure distribution prevailing at the nth 
iteration is used in the solution of the axial and radial velocity momentuc 
equations which can be written in the general form: 


m m 


im 


(4.4.1) 


where U. represents either of the velocity components, u or v, and the 
operator is a difference operator given by 


A. 

ip 


P„ ■ Pe 


for * u 


■ P* ’ Pn 


for 


(4.4.2) 
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The subscript m denotes the nodes, E, W, NE, 

Fi^s^ Corrector Sta£^ - A new velocity field u, ** tosether .> 
r«po„<l.ng „ev pr„.u„ field p. clcuLt.d .u?h ‘ 
divergence (incompressible continuity) condition 


du. 

1 


( 4 . 4 . 3 ) 


is satisfied. The momentum equation is now written as: 

“i “ m “im ‘ (4.4.4) 

i. "’■ • v.locit,-ipcre.,„. 

V • ★* ^ ... a. 

A ^ (u. “ u.) « D ^ Z) . (p / 

P 1 1 p '■i 'P P (4.4.5) 

Combining the divergence of this equation with Eq. (4.4.3) the 
pressure-correction equation is obtained: 


i l^p' ^ ^p' ]p' ' “i 


( 4 . 4.6 


where 


p* - p" 


( 4 . 4 . 7 ) 


revised pressure distri- 

tion p 16 determined and a new velocity field u^** is calculated. 

Secoad_Corrector_Sie£ - A new velocity field u-*** together with a 
corresponding pressure distribution p** are calculated such that 


diUi - 0 

The momentum equation is now written as: 

n • ** 


( 4 . 4 . 6 ; 


.u; 


u . 


r A 

m in 


‘ - “p* 


♦ s'"! 


P 1 

Subtracting Eq. (4.4.5) from Eq. (4.4.9) yields: 

*p '“i ^ *:■ (-i -i). - bi(p**-p*) 


( 4 . 4 . 9 ) 


(4.4.10) 


“ I I, I I 


(4.4.11 ) 


where 
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(A. 4. 12) 


Eq. (4.4.11) is solved for the p” field so that the revised pressure dis- 
tribution p** is computed and a new velocity field u.*** is determined. 

The pressure and velocity distributions for the (n+17th iteration are taken 
as p** and respectively. 

A similar procedure is used for the turbulent kinetic energy and 
turbulent energy dissipation equations because these equations are strongly 
coupled by means of their source terms. The equations for k and e are 
first solved for the source term components based on the previous iteration 
values. A corrector step is then applied to account for the changes in 
the source terms. 
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5.0 DESCRIPTION OF COMPUTER PROGRAM 


The computer program based on the mathematical model described in 
Section 3.0 is described in this section. Information is presented to 
assist the user in the installation of the program on modern large-scale 
computers. Additionally, an input format is provided. 


5.1 General Remarks 


The computer program is intended for use in the batch-operating mode 
on modern large-scale computers. Specifically, programming constructions 
have been avoided that refer to the characteristics of specific computer 
operating systems. In addition, all logical unite referred to by the 
program are defined using integer constants initialized in the BLOCK DATA 
routine BLOCKD. In the present version of the program, these logical units 
have been assigned as follows: 


Name 

Value 

Puroose 

LU4 

4 

Restart file 

LU5 

5 

Card reader 

LU6 

6 

Line printer 

LU6 

8 

Scratch unit 

LU9 

9 

Scratch unit 


The progr^ may be stopped periodically and the results stored in the 
restart file, LDA. The calculation may subsequently be continued from the 
latest Iteration by using the restart file in conjunction with the options 
described on Card No. 2 of the Input Format, Section 5.7. 

The present program was derived from the development version provided 
College to UTRC. The Imperial College computer is taken to be 
a CDC 6600 computer with a FORTRAN IV compiler. Except for eliminating tvo 
minor incompatible constructions and the CDC PROGRAM card, the develop- 
ment program was completely compatible with the FORTRAN V compiler used on 
the UNIVAC 1100/81A computer at UTRC. For the convenience of users of CDC 
computers, the PROGRAM card has been retained but rendered inoperable in 
the present version (see routine AMAIN, lines 1-2). It is believed that 
the present version is also compatible with the IBM FORTRAN IV compiler. 

No attempt was made to compile the present version of the program 
using the UNIVAC 1977 ASCII FORTRAN compiler since it would have been 
necessary to eliminate any incompatibilities between the constructions used 
in the development version and the 1977 ASCII standard; in addition, it is 

known that some 1977 ASCII standard compilers are not compatible with each 
other. 

The workers at Imperial College have indicated throughout the program 
those sections of the code that are restricted to flow within a research 
combustor of the type shown in Fig. 1. They have labeled these sections of 
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the code as "Problem Dependeat"* No attempt has been made to determine if 
other sections of the code should be designated in a similar manner. 

The computer program can be allocated in approximately 50000 decimal 
words, or 200000 bytes. 


5.2 Program Flow Chart 

A general flow chart of the computer program is shown in Fig. 6. The 
input to the program is read in subroutine MODINP; additional constants are 
initialized in the BLOCK DATA routine. Initial guesses are assigned for 
the distributions of each of the field variables. The program then 
proceeds as follows: 


1. A pressure distribution is guessed. For the first iteration, 
the estimated pressure distribution is provided by subroutine 
MODINP. On subsequent iterations, it is assumed to be the 
distribution determined at the end of the previous iteration. 

2. The axial and radial momentum equations are solved for the 
axial and radial velocity distributions using subroutines 
COEFU and COEFV. 

3. The pressure-correction equations (i.e., the PISO algorithm) 
are solved using subroutine CALUVP to obtain the corrected 
pressure, axial velcoity and radial velocity distributions. 

4. The swirl velocity (actually, the product of radius and swirl 
velocity) distribution is solved for using subroutine CALCFI. 

5. The turbulent kinetic energy and energy dissipation rate are 
determined using subroutine CALCFI. 

6. The turbulent kinetic energy and energy dissipation rate are 
adjusted using the predictor-corrector technique noted in 
Section 4.4. 

7. The distributions of the remaining dependent variables (e.g., 
enthalpy, mixture fraction, fuel concentration) are 
calculated using subroutine CALCFI. 

8. The absolute values of the normalized residuals for each 
variable are examined at all points in the flow field. If 
the residuals are all less than a apecified maximum value, 
the solution is considered to be converged; otherwise, the 
iteration is repeated from Step 1. 


The distributions for each variable are solved using a line-by-line itera- 
tive procedure which utilizes a tridiagonal matrix algorithm (subroutine 
LISOLV) to solve simultaneously for the variables along each grid line. An 
alternating direction scanning procedure is employed to sweep the whole 
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field. The number of eveeps is usually different for each variable. Gene- 
rally, two aveeps are used for the axial and radial velocities and one 
sweep is used for all other variables except for the primary pressure 
correction equation. Since this equation is strongly elliptic and has 
Neumann boundary conditions specified at all boundaries, at least three 
sweeps are generally made. 


3.3 Description of Routines 

A brief description of each of the routines used in the computer pro- 
gram is presented in this section. The principal FORTRAN variables are 
listed in Section 5.4, the entry point and external reference table is 
presented in Section 3.3, the source language listings are displayed in 
Section 5.6, and the input format is presented in Section 5.7. 

AMAIN 


AMAIN is the main program whose principal operation is to control the 
progress of the iterative solution described in the previous section. 
Additionally, AMAIN controls the rate of intermediate printout (e.g., the 
printing of the maximum residuals and values of the dependent variables at 
the monitoring grid node) and it determines which of the distributions for 
the dependent variables will be displayed by subroutine PRINT. For react- 
ing flows, this routine initiates the combustion reaction by introducing at 
the iteration ITREAX an artifical difference in the mixture fraction and 
fuel concentration (e.g., by setting m^ ■ 0.9f) at the monitoring grid node 
(IMON, JMON). Finally, AMAIN retrieves dependent variable distributions 
from the file on logical unit LU4 to restart a calculation; and it stores 
these distributions in the same file at the end of the current execution so 
that the case may be started subsequently from the latest iteration. 

BLOCKD 

The BLOCK DATA routine initializes certain print heading arrays, 
assigns values to the constants used for logical units and for the 
turbulence model, and sets some miscellaneous controls. 

BOUNDS 


This subroutine determines the values of the flux blending factor 
used to adjust the flux contributions to the coefficients of the finite- 
difference equations such that the coefficients are non— negative. The 
bounding procedure is described in Section 4.1.4. It should be noted that 
Y does not appear explicitly in this subroutine; instead, the values of 
the flux contributions are reduced in accordance with the hounding proce- 
dure. 

CALCFI 


Subroutine CALCFI is a general subroutine for calculating the distri- 
butions of all scalar variables and the product of swirl velocity with 
radial position. The coefficients of the finite— difference equations are 
calculated for each variable (i.e., for each call to CALCFI). In some 
cases, it is known that the co— efficients for one scalar are simply related 
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to those for Another* A provision is ude to determine one set of co- 
efficients from the other without the need to recompute all of the flux 
contributions. Often, two sets have identical values or differ only by the 
ratio of the Frandtl numbers for each variable. The source term contribu- 
tions to the finite-difference equations are then calculated by calling 
subroutine SORCFI. The modifications to the co-efficients due to the 
boundary conditions are determined by calling subroutine MODFl. The co- 
efficients are then bounded using subroutine BOUNDS. The underrelaxation 
factors are then applied. The system of simultaneous equations is 
then solved using subroutine LISOLV. 

CALUVP 


Subroutine CALUVF uses the coefficients assembled in subroutines COEFU 
and COEFV for the axial and radial momentum equations to solve for the 
initial estimates of the axial and radial velocity distributions, 
respectively, for the current iteration. The co-efficients for the finite- 
difference form of the primary and secondary pressure-correction 
(predictor— corrector ) equations are then calculated and the final velocity 
fields are determined. (See the discussion of the FISO algorithm. Section 
4 . 4 .) 

COEFU 


Subroutine COEFU is used to calculate the co-efficients of the finite- 
difference form of the axial momentum equation. First, the convection and 
diffusion contributions to the coefficients are calculated. Next, the com- 
ponents of the source term are determined. The flux contributions to the 
finite-difference coefficients are then calculated. The effects of the 
boundary conditions are incorporated by calling the entry point MODU of 
subroutine MODUVF. The bounding strategy is then applied by using 
subroutine BOUNDS. Finally, the coefficients are stored for use by sub- 
routine CALUVP. 

COEFV 


Subroutine COEFV is used for calculating the coefficients for the 
finite-difference form of the radial momentum equation. Its operation is 
similar to that of subroutine COEFU. The boundary conditions are applied 
by calling the entry point MODV of subroutine MODUVF. 

CORECT 


This subroutine determines the values of the flux blending factor used 
by subroutine BOUNDS to limit the co-efficients of the finite-difference 
equations to non-negative values. 

PORTED 


The turbulent kinetic energy and energy dissipation rate equations are 
strongly coupled to each other by their source terms. Subroutine PORTED 
performs a predictor-corrector operation for the two equations to provide 
more accurate estimates of these variables and thereby improve the rate of 
convergence of the entire system of equations. 
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INIT 


Subroutine INIT defines the grids used in the calculation, determines 
cell dimensions for the control volumes, and calculates the factors used 
for interpolating variables at each face of the control volumes. The sub- 
routine also initializes all arrays to appropriate values. 

LISOLV 


This subroutine solves the system of simultaneous finite-difference 
equations for each variable using an alternating direction method with a 
tridiagonal matrix solution algorithm. The first sweep of the system of 
equations is from vest to east, the second form south to north, etc. The 
total number of sweeps is specified for each variable in the input. 

MODFI 


Subroutine MODFI is used to modify the coefficients of the finite- 
difference equations for each scalar variable due to the influence of the 
boundary conditions appropriate for that variable. The scalars are: tur- 
bulent kinetic energy, turbulent energy dissipation rate, enthalpy, fuel 
mass fraction, mixture fraction, and swirl velocity (actually, the product 
of swirl velocity with radial position). The boundary conditions are 
described in Section A. 3. 

MODIKP 


This subroutine reads the input cards, writes the input on the line 
printer, sets all logical variables, calculates the reference values for 
the residuals, and initializes the dependent variables. 

MODUVP 

Subroutine MODUVP modifies the coefficients for the axial and radial 
momentum equations and for the pressure-correction coefficients in accor- 
dance with the boundary conditions (Section 4.3). The entry point MODPRO 
was provided by the authors of the program at Imperial College to deal with 
the effects of transport properties on the coefficients that are not incor- 
porated elsewhere in the program; presently, this entry point contains no 
executable code. 

OUTPUT 


Subroutine OUTPUT prints the values of selected integrated quantities 
(main flow rate, secondary flow rate, etc.) and other parameters that 
characterize the flow field. For reacting flows, the distributions of 
oxygen, nitrogen, water vapor, and carbon dioxide are also printed. 

PRINT 


This subroutine generates a printout of the tabulated values of the 
field variables. 
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PROPS 


Subroutine PROPS is used to determine the fluid properties at each 
point in the flow field. It is assumed that the main flow contains air 
with an oxygen mass fraction of 0.233 and a nitrogen mass fraction of 
0.767. The average heat capacity, molecular weight, and density can then 
be determined. The viscosity at each point in the flow field is the sum 
of the laminar viscosity and turbulent viscosity. The effective 
diffusivity for the turbulent kinetic energy equation is also calculated. 
The diffusivities for other scalars are determined in subroutine CALCFI 
by multiplying the diffusivity for the turbulent kinetic energy equation by 
the ratio of the turbulent Prandtl number for another scalar to that for 
the turbulent kinetic energy. 

SORCFI 


This subroutine calculates the linearized source term components for 
each scalar in accordance with Section 4.2. 


5.4 Principal FORTRAN Variables 


The principal FORTRAN variables are listed in this section. The 
variables are arranged in logically occurring groups rather than in strict 
alphabetical order. Variables written as var(I,J) denote two-dimensional 
arrays defined for each axial gridline location I and each radial gridline 
location J. 


AE(I,J) 

AW(I,J) 

AN(I,J) 

AS(I,J) 

ANE(I,J) 

ANW(1,J) 

ASE(I,J) 

ASW(1,J) 

AP(I.J) 

CAPPA,EL0G 

CCEN,CIN 


CMU,CD,C1,C2 

CPM(I,J) 

CP02 

CPN2 

CPC02 

CPH20 

CPF 


Coefficients for the finite-difference equation for the 
control volume at point (I,J). 

E ■ east, W - west, etc., P - point (I,J) 


Constants in the velocity law of the wall 

Input mass fraction of second species (fuel) in secondary 
and main stream, respectively. 

Constants in the two-equation turbulence model 

Specific heat of mixture 

Specific heats of oxygen, nitrogen, carbon dioxide, water 
vapor and fuel, respectively 


DELCEN, DELTA 


Input value of boundary layer thicknesses for secondary 
and main stream, respectively 
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DEN(I,J) 

DENS IT 

ED(I,J) 

EN(I,J) 

EKI.J) 

E2(I,J) 

E3(1,J) 

F(1,J) 

FM(I,J) 

FN1(1,J) 

FN2(I,J) 

FN3(I,J) 

F02 

FN2 

FC02 

FH20 

FPR 

GAMH(I.J) 

HR 

lED 

lEN 

IF 

IFM 

IP 

IPP 

ISWR 

ITE 

lU 

IV 

IMON.JMON 

INCAL 


INCOMP 


INDCOS 


Density 

Input value for density for incompressible flow case 
Turbulent kinetic energy dissipation rate 
Enthalpy 

Flux contributions, east face of control volume 

Mixture fraction 
Fuel mass fraction 

Flux contributions, north face of control volume 


Mass fraction of oxygen, nitrogen, carbon dioxide, water 
vapor, and products of reaction, respectively 


Diffusivity for turbulent kinetic energy 
Heat of reaction 

Parameters used to set the logical array INCAL for turbu- 
lent energy dissipation rate, enthalpy, mixture fraction, 
fuel mass fraction, primary and secondary pressure 
corrections, swirl velocity, turbulent kinetic energy, 
axial velocity, and radial velocity, respectively 


Location for grid node used to monitor the calculation 

Linear array for the logical variable that determines if 
the differential equation for the Kth dependent variable 
is to be solved (1NCAL(R)>.TRDE. ) or not solved 
(INCAL(K)«. FALSE.) 

Logical variable indicating whether flow is compressible 
(.FALSE.) or incompressible (.TRUE.) 

Parameter to indicate whether flow is two-dimensional (1) 
or axisymmetric (2) 
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INDPRI 

INPRO 

IPREF.JPREF 

ISCEME 

IT,JT 

ITREAX 

ITSTEP 

IVISCO 

JEXIT 

JINS 

MAXIT 

MAXSVP 

MODOP 

NI.NJ 

NSUP 

NUMPRI 

P(I.J) 

PERR 

PIN 

PRANDL 


The nmiiber of iterations between the output of the tabu- 
lated field variables 

Logical variable indicating whether flow properties are 
constant (.FALSE.) or variable (.TRUE.) 

Location of reference pressure node 

Parameter indicating whether hybrid differencing (1) or 
skew-upwind differencing (2) is to be used 

Maximum dimensions of the arrays var (I,J) 

Iteration at which chemical reaction starts 

Number of iterations to be performed for current 
execution 

Logical variable indicating whether viscosity is constant 
(.IRUE.) or variable (.FALSE.) 

Radial grid line number for location of maximum radius 
at exit 

Radial grid line number for location of maximum radius 
for secondary inlet 

The maximum number of iterations to be executed 

The maximum number of sweeps to be executed by LISOLV for 
solvii^ the pressure-correction equation 

Parameter indicating whether a new case is being started 
(1) or a previous case is being run using the Restart 
file (2) 

Number of input axial gridline and radial grid line 
locations, respectively 

A linear array denoting the number of sweeps to be used 
when solving for each of the dependent variables 

The nnmber of iterations between output of the maximum 
residuals and monitoring node information 

Pressure 

Maximum allowable residual in pressure correction 
equation 

Input vxlue of pressure 
Laminar Prandtl number 
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PRANDT 

PREF 

PP(1,J) 

SORHAX 

SHORN 

SP(I.J) 

SU(I,J) 

STOIC 

T(I,J) 

TE(l.J) 

TCEK.TIN 

TURBCN, TURBIN 

U(1,J) 

UCEN.UIK 

UGC 

URF 

URFDEN 

URFVIS 

V(I.J) 

VIS(I.J) 

VISCOS 

W(I,J) 

WCEN.WIN 


A linear array of turbulent Prandtl numbers for each 
dependent variable 

Reference pressure at node ( IPREF , JPREF ) 

Pressure correction 

Maximum allowable residual for any dependent variable 

A linear array of the maximum residuals for the dependent 
variables 

Coefficients in linearized source term that are dependent 
on and independent of the dependent variable at (I,J), 
respectively 

Stoichiometric oxygen to fuel mass ratio 
Temperature 

Turbulent kinetic energy 

Input values of temperature for secondary and mainstream, 
respectively 

Input values of constant of proportionality for turbulent 
kinetic energy for secondary and main stream, 
respectively 

Axial velocity 

Input values of axial velocity for secondary and main 
stream, respectively 

Universal gas constant 

Linear array for underrelaxation factors for dependent 
variables except for density and viscosity 

Underrelaxation factor for density 

Underrelazation factor for viscosity 

Radial velocity 

Total (laminar plus turbulent) viscosity 
Input value of laminar viscosity 
Swirl velocity 

Input values of swirl to axial velocity ratio for 
secondary and main stream, respectively 
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XU Input linear array of axial locations of axial velocity 

grid lines 

YV Input linear array of radial locations of radial velocity 

grid lines 

5.5 Entry Points and External References 

The entry points for each routine and the external references by each 

routine are listed in this section. An asterisk is used to denote external 

references to FORTRAN library routines. 

Routine Entry Points External References 

AMAIN AMAIN AMAXl* 

CALCFI 

CALUVP 

COEFU 

CALFV 

CORTED 

GAM 

INITOP 

INTOP 

PARAM 

PRINT 

PROPS 

STARTV 

BOUNDS BOUNDS CORECT 

CALCFI CALCFI AMINI* 

BOUNDS 

LISOLV 

MODFI 

SIGN* 

SORCFI 

CALUVP CALUVP LISOLV 

MODP 

MODUl 

COEFU COEFU AMAXl* 

AMINI* 

BOUNDS 

MODU 

SIGN* 

COEFV COEFV AMAXl* 

AMINI* 

BOUNDS 

MODV 

SIGN* 

CORECT CORECT AMAXl* 

AMINI* 
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CORTED 


CORTED 


MODFI 

SORCFl 

Vise 


INIT grid 

INIT 

SET 


LISOLV 

LISOLV 

MODFI 

MODFI 

MODIKP 

MODINP 


PARAM 


STARTV 

MODUVP 

MODP 


MODPRO 


MODU 


MODUVP 


MODUl 


MODV 

OUTPUT 

FINOP 


INITOP 


INTOP 


OUTPUT 

PRINT 

PRINT 

PROPS 

DENS 


GAK 


PROPS 


Vise 

SORCFl 

SORCFl 


AMAXl* 

GRID 

PROPS 

SET 


PRINT 

MINO* 

AMAXl* 

AMAXl* 

AMINl* 
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5.6 Source Listings of Routines 


The listing of the FORTRAN source language for each routine is presen- 
ted in this section. 
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ITREAX Start chemical reaction on iteration ITREAX; this variable is used to permit flow calculation to 

stabilize prior to start of combustion. Typically, ITREAX = 30. 




6.0 TEST CASES 


Three test cases are described in this aection. The first two cases 
are based upon flow conditions and data supplied to UTRC by W. M. Roque- 
of Wright-Pat terson Air Force Base in consultation with the NASA 
Monitor. These data were obtained for nonswirling flows in a 
research combustor of the type shown in Fig. 1. Both a cold flow run and a 
hot flow run are simulated. The cold flow and hot flow nonswirling test 
cases are designated cases WPAFBl and WPAFB2, respectively. The third case 
(WPAFB3) xs derived from the hot flow test case without swirl by including 
awxrl in the main (outer) stream and neglecting chemical reaction. It is 
included only to provide confirmation that the awirl option is operational. 

•11 the cases, the main flow consisted of air with a flow rate of 
2 kg/sec. For the cold flow test cases, the secondary flow consisted of 
carbon dioxide having a flowrate of 6 kg/hr; for the hot flow test case, 
the same flow rate was used but the carbon dioxide was replaced with 
propane. The secondary inlet temperature and flow velocity also differed 
between the hot and cold flow cases. It was assumed that the turbulence 
was isotropic. From the test data provided to UTRC, it was estimated that 
the ratio of turbulent kinetic energy to the square of the flow velocity in 

each stream was 0.003. The inlet flow conditions are summarised in Table 
7 • 


The primary purpose of the test cases was to demonstrate that the 
computer program is capable of producing a physically realistic result w’ith 
a minimum of user intervention. Achieving an accurate comparison between 
calculated and measured results was considered to be of secondary impor- 
tance. Therefore, a relatively coarse grid system, consisting of twelve 
axial and ten radial grid lines, was used for all calculations. The grid 
system is shown schematically in Fig. 7. The numbering of the grid lines 
xs consistent with that defined in Sections 5.4 and 5.7. 
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6.1 Cold Flov Test Cate Without Swirl 

The cold flow test cete without twirl (ctte WPAFBl) wet run for t 
total of 100 iterationt. A fitting of the input cardt and a copy of the 
printout corretponding to iterationt 1 through 5 “ 

Appendix 1. A copy of the printout for iterationt 31 through 100 it alto 

included in AppendU 1. All input waluet for underrelaxatiw 

number of tweept per iteration, etc, were the tame for all of the runt. 

At the end of 100 iterationt, the calculation had not converged tince 

tome of the retidualt exceeded the input “ooS^^^The retid- 

Cenerally, SOEMAX it attigned a value between 0.001 

ualt for fuel concentration (F.CONC.) and mixture fraction (M.FRAC.) are 
about 0.01. However, for purpotet of illuatrating the operation of the 
program, thete retidualt may be contidered to be tmall Mough. Ho •“^Jable 
length tcale hat been defined in the program for normalizing the turbule t 
energy dittipation rate. It it teen that the final value of the retidual 
(DISS) for thit variable it 1.78; however, tince the value of DISS on the 
firtt iteration wat 96200, the iteration for turbulent raergy dittipation 
hat probably converged. Therefore, the cate hat 

of the 100th iteration. Approximately 1670 aecondt of UNIVAC 1100/81A 
machine time were uted for thit calculation. 

The computed variation of centerline velocity with axial dittance is 
thown by the tolid line in Fig. 8. The meatured data corretpond to carbon 
dioxide flow ratet of 4 and 8 kg/hr; data are not available for the carbon 
dioxide flow rate of 6 kg/hr uaed in the calculation. It can be teen that 
the calculated reaultt show a more rapid initial deceleration of centerline 
velocity than that ahown by the data. Thit could be the retult of a 
greater rate of mixing of the two streams with the rapid mixing the result 
of severe numerical diffusion due to the very coarse grid used for the 
calculation. Similar behavior is exhibited by the calculated results for 
axial variation of centerline carbon dioxide concentration thown in Fig. V 
which are compared to measured data at the simulated flow rate. 
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6.2 Hot Flow Test Case Without Swirl 

without swirl (case WPAFB2) was run for 125 
Iterations, a listing of the input cards and a copy of the results for 
Iteration. 1 throngh 125 .ro .ho™ in Appendix 2. %he iiLrje.ILfriot 
the the mixture fraction is slightly greater than 

«»intiL“^.r'T'” ‘ *<>' ‘X‘ torbulenr 

. «Urof 463 on >35000 on the fir.t iteration tc 

Use ofth^ ^63 on the last iteration. Therefore, this case has converged 
e of the CMbustion model was initiated on the 30th iteration All ^ 

”«rtJ™;™rt';r°”-”“""" »» ^—uon! 

aecond. of OKIVAC IlOO/iu'naculle'u^^''’"'' *'"’'‘>“““>3 3080 

Fie in^ calculated and measured variations of axial velocity are shown ir 
Fig. 10 and of temperature in Fig. 11. a comparison of calculated and 

rh! predicted results exhibit greater mixing between 

the ^in and secondary streams than that shown by the data- again the 

^errcLrsrgrid"'' numerical diffusion caused by the use of a 
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6.3 Cold Flow Test Case With Swirl 

A cold flow test case (case WPAFB3) was run to demonstrate that the 
computer program was operating properly. The case was terminated after 
twenty iterations. A listing of the input cards and a copy of the results 
are included in Appendix 3. The calculation has obviously not converged 
but appears to be stable. This test case is provided only for documenta- 
tion purposes and to aid the user in installing the computer program on 
another computer system. 
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8.0 LIST OF SYMBOLS 


A 

a 

a 

C 

S 

,c^ 

D 


E1.E2.E3 

F 

f 

G 

h 

Hr 

1 

i 

J 

K 

k 

M 

n 

N1,N2,N3 

Pe 


Rf 

S 


Co~efficient in f inite~difference equations 

Flow area 

Flux contribution 

Convective flux 

Heat capacity 

Constants in turbulence model 

Diffusive flux 

Constant in lav of the wall 

Flux contributions, east face of control volume 
Flux 

Mixture fraction 
Turbulence generation tenn 
Enthalpy 

Heat of reaction 

Index for nodes - axial direction 
Stoichiometric oxygen to fuel mass ratio 
Index for nodes — radial direction 
Constant in law of the wall 

Skewing interpolation factor or turbulent kinetic energy 
Molecular weight 
Mass fraction 

Flux contributions, north face of control volume 

Pec let number 

Universal gas constant 

Reaction rate 

Source 
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s 

u 

V 

X 

y 

3 


Y 

£ 


:v 


cx.iy 

lx.ly 


p 


,u 


c 


(t 


Distance along streamline 
Axial velocity 
Radial velocity 
Axial co-ordinate 
Distance from vail 
Interpolation factor 
Exchange co-efficient 
Blending factor 

Coordinate indicator (0 for planar, 1 for axisymmetric) 

Differential volume 

Control volume dimensions 

Internodal distances 

Turbulent energy dissipation rate 

Hybrid differencing switch (1 for central differencing, 0 for 
bounded skew-upwind differencing) 

Dependent variable 

Density 

Flow velocity direction switches 
Turbulent Prandtl number for (t 
Shear stress 


Subscripts 

BSHD 

BSUD 

CD 

CO2 

E 

e 

eff 


Bounded skew hybrid differencing 

Bounded skew-upwind differencing 

Central differencing 

Carbon dioxide 

East node 

East face 

Effective value 
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f 

H20 


ID 

max 

min 

N 

K2 

n 

KE 

NV 

P 


pr 

5 

6 

SE 

sw 

SUD 

t 

u 

UD 

W 

V 

< 


Fuel 

Water vapor 
Laminar 
Mean value 
Maximum value 
Minimum value 
North node 
Nitrogen 
North face 
Northeast node 
Northwest node 
Oxygen 

Node at center of control volume 

Product 

South node 

South face 

Southeast node 

Southwest node 

Skew upwind differencing 

Turbulent 

Value not dependent on 
Upwind differencing 
Vest node 

West face or twirl velocity 
Dependent variable 


Superscripts 

( 


Axial velocity 
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Radial velocity 


( )* 
( )• 
( )” 


Transition 

For co-efficients, value not associated with 
For co-efficients, value associated with 
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TABLE 1 


Turbulence Model Constants 


Constant 


Value 


Cl 




E 

K 


1.44 
1.92 
0.09 
1 .0 

k 2/|(C2 - Cl) 
9.0 

0.4187 
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TABLE 2 

Constants in the Polynomials for Heat 
Capacity for Various Species 

^2° C3H8 O2 CCp 

^^1 2.7168 -0.24515 3.620 3.8963 \ AiZl 

C 2 2.9454>10‘3 0.36714x10"1 7.3612x10"^ 1.5155xl0‘*= 3.0?£2>i:’' 

C 3 -8.0224x10-7 -0.19576x10*4 -1.9652x10*7 -5.7235x10-7 -1.2392>1C-^ 

C 4 1.0227X10-10 0.46527x10-2 3.6202xl0-''l 9. 9807x1 0-l ^ 2 . , v- -I : 

C 5 -4.8472x10-"'2 -0.29008x10-''2 -2.8946x1 O-^ 2 -6.5224x10'''^ -..552c -I'-' 

2l 4.0701 -0.24515 3.6256 3.6748 2.£0:i 

C 2 -1.1084x10-3 0.36714x10-'' -1.8782x10-3 -1.2081x10-3 8.7351 >tc-3 

C 3 4.1521x10-2 -0.19576x'0-4 7.0555x10-^ 2.324x10-2 

I C 4 -2.9637x10-- 0.46527x10-2 -6.7635xl0’2 -6.3218x10-^2 2 . 0.:22 > ' : ‘ ' 

C 5 8.0702x10-''3 ~0.2900Bx10-''3 2.1556x10-^2 -2 . 2577x1 n-'' 3 


9-3 







*1 


TABLE 3 

Definition of Bounded Skew Hybrid Differencing Fluxes 


Vest 


F 


'-’BSHD 




V w, 


CD 


(1 ~ > ) F 
w 


'^BSLT) 


East 


^BSHD 




(1 


) T 


South 


DS 


's^scD 


+ (1 


' ) F 

"^SBSUD 


North 


Dn 




+ (1 


> 1 F 

"bSUD 
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TABLE A 


Central Differencing Fluxes 


We 5 1 




‘CP 


Dv 


Pe.. (1 - a„) ^ 1 


V. 




East 


"CD 


Pe„ (1 - 4 1 


C 

P 


^-^e^ee “ ^ ) 


South 


CD 


Pe. n 


'» 4 1 


( ' sPpj ” 0 


North 


"CD 


Pe_ (1 - aj 4 1 


-n n 


1 'P * ( ■ ‘ ) ‘n 
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TABLE 5 


Bounded Skew-Upwind Differencing Fluxes 


Vest 


U : 1 ^.:-^ = sw -»-<> ♦'^’1 


(1 - c“) , -p 

V I y 


‘ V ^ S 


, V, I 

- » - "v> ‘N 


East 


'eBSUD 


= ?e [ -P 
e y 


(1 - r^) :r’ ' ^Pe„ ' Pe,> •' 


_u 

e 


'P " e 'S 


(1 


-O =J 


+ (1 - cp 

e 


- % 'SE - <1 - 


e KT 


South 


■®BSIT) _ ^ r V , _ _Vx , 1 _ (p - Fg ) Y k, 

= Peg Kg "g ^ S'* -p J ®s eg s £ 


s I 's " s 


- -""w - c“ d I 

'SE s SV I 


(1 - j ^ 


-% Y .|1 


North 


"BSUD 

I>N 




% 1 'p ■ “ " p ’ '" K '- ' 


♦ (1 - cj - n 


“11 

-c c ^1 

n NI n NVj J 
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TABLE 6 


Co-e f f i c ien t s of the Finite-Dif ference Equation? 


A^.(I,J) = E2(I-1,J) - E3)j_^^j 

+ >s [(1 - r'") N3]^ , . 

^ n n •'I.J-l 

- ^ri [ r c N3U , 
n n i J 


A^d.J) = E1(I,J) + [(1 - E3 ] , + [(1 - c'^Kl - K31 

® ^ ® n n ' I , 


(1 - K 3 ^ 

nr. r. ^1,J 


Ac(i,J> = K2(I .J-1) - 


■ N3] 4 -> f(i 


r'') c'* E 3 ^ 
e e I“1,J 


u \* , 

' r E3 
€ e I .j 


Aj;(l,J) = NKI.J) 4 [(1 - c') K3] 4 v^, [(1 - c"^)(l - r'') E3] 

n i , j e e I- ] 


e (1 - 13] 

^ e e ^l.J 




Ag£ ( I , J ) 

Aj^r;(l , J) 
Ake(I.J) 


: (1 - c ) K3] - 1(1 - c“) c E3)^ ^ 

n ri i , J-i e e I , J 


V '-e - %> c“ K3l^_, 


[(1 - r )(1 - c'’) E3] 


e 'e''" 'e' ""'l J ” - c )(1 - -“)]N3 

^ ^ ^ n n I,J 


NOTE: The blending factors Y^ and y„ are the factors appropr 

for the (I,J) node as determined by the procedure discussed in 
Section 4.1.4. 
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TABLE 7 


Inlet 

Flow Conditions 

for Test Cases 

Case 

WPAFBl 

WPAFB2 

Ma in_F_low 



Axial Velocity, tn/sec 

49.3 

49.3 

Swirl/Axial Velocity 

0.0 

0.0 

Tenpera tur e , K 

293 

293 

Turbulent Kinetic Energy/ 

0.005 

0.005 

(Main Velocity)^ 



Secondarv Flow 



Axial Velocity, m/sec 

52.5 

69.6 

Svrirl/Axial Velocity 

0.0 

0.0 

Temperature, K 

293 

400 

Turbulent Kinetij; Energy/ 

0.005 

0.005 

(Sec. Velocity)" 



Pressure, Pa 

98000 

98000 


VsTAFBS 

49.3 

0.101 

293 

0.005 

69.6 

0.0 

400 

0.005 

96000 
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Fig. 1 Geometry of Bluff-Body Research Combustor 





Fig. 2 Grid System for Scalars and Tangential Velocity 
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Fig. 4 Control Volume for Scalar <h (I, J) 
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note indices I, J REFER TO LOCATION OF SCALAR 

Fig. 5 Control Volume for Skew Upwind Differencing for West Face of Control Volume 
and Positive Velocity Components 
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S3-8-30-12 
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Fig. 7 Grid System Used for Test Cases 



COj flowrate = 6 Kg/rif 



Fig. 9 Centtrlln® CO 2 Concentration Distribution for Cold Flow Test Case Without Swirl 
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TEMPERATURE — deg K 


Fig. 12 Temperature Profile for Hot Flow Test Case Without Swirl 
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